Variational Cluster Perturbation Theory for Bose-Hubbard models 
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We discuss the application of the variational cluster perturbation theory (VCPT) to the Mott- 
insulator-to-superfluid transition in the Bose-Hubbard model. We show how the VCPT can be 
formulated in such a way that it gives a translation invariant excitation spectrum - free of spuri- 
ous gaps - despite the fact that if formally breaks translation invariance. The phase diagram and 
the single-particle Green function in the insulating phase are obtained for one-dimensional systems. 
When the chemical potential of the cluster is taken as a variational parameter, the VCPT reproduces 
the dimension dependence of the phase diagram even for one-site clusters. We find a good quan- 
titative agreement with the results of the density-matrix renormalization group when the number 
of sites in the cluster becomes of order 10. The extension of the method to the superfluid phase is 
discussed. 

PACS numbers: 05.30.Jp, 73.43. Nq, 03.75.Lm 



I, INTRODUCTION 

The Bose-Hubbard model describes interacting bosons 
on a lattice. It provides a generic description of the 
quantum phase transition between superfluid (SF) and 
Mott-insulator (MI) states observed in condensed-matter 
systems such as Josephson junction arrays or granular 
superconductors^ as well as in ultracold atoms in opti- 
cal latticesi 2 ^*^ The remarkable degree of experimental 
control over all the relevant parameters (density, inter- 
action strength, lattice geometry and dimensionality) in 
ultracold atoms makes possible a detailed study of the 
MI-SF transition. 

The Bose-Hubbard model has been studied nu- 
merically using the Gutzwiller mean-field ansatzpii^ 
the density-matrix renormalization group(2*ifi exact 
diagonalizations ) 11 ' 12 and quantum Monte Carloi 1 ^ 14 ' 1 !? 
(More recent works include the harmonic trap poten- 
tial that confines the ultracold atomic gases; see for in- 
stance Ref. 0.) Most analytical approaches rely on a 
perturbation theory that assumes the kinetic energy to 
be small and treats exactly the on-site repulsion. The 
intersite hopping is taken into account either at the 
mean-field level or in perturbation in a strong-coupling 
expansion, 1 - 1 7, 1 V P, 2 9^i ,22,23,24,25,26,27 

In this paper, we apply the variational cluster pertur- 
bation theory (VCPT) to the Bose-Hubbard model at 
commensurate density. The VCPT has been developed 
for strongly-correlated fermion systems. It is an exten- 
sion of cluster perturbation theory (CPT) that is based 
on the self-energy functional approach (SFA). Within 
the CPT) 28 i 29 i 30 the lattice is partitioned into discon- 
nected identical clusters. The Hamiltonian of the clus- 
ter is solved numerically, and the intercluster hopping is 
then treated perturbatively to leading order in a strong- 
coupling expansion. Contrary to exact diagonalizations 
of small systems, the CPT provides results in the thermo- 
dynamic limit, and the single-particle Green function is 
defined for any wave vector in the Brillouin zone. When 



based on a single-site cluster, the CPT yields the Hub- 
bard I approximation j2l applied to bosonic systems, it re- 
produces the leading order of the aforementioned strong- 
coupling theory^ 

The SFA is based on the variational principle 
<5Q[£]/5£ = for the grand potential expressed as a 
functional of the self-energy Y]M£2*& The stationary con- 
dition is solved within a restricted space of self-energies 
taken from a reference system that can be solved nu- 
merically. Within the VCPT^^ the reference system 
consists of disconnected identical clusters. The VCPT 
improves on the CPT since the parameters of the intr- 
acluster kinetic Hamiltonian are variational. In partic- 
ular, this enables to consider broken-symmetry states. 
The VCPT has been used with some success to study 
strongly correlated electron systems ,^^37,^39 

The outline of the paper is as follows. In Sec. |n] we 
describe the SFA and the VCPT for bosonic models in 
the absence of superfluidity. We modify the original for- 
mulation of the SFA2 2 " to ensure that final results remain 
translation invariant regardless of the choice of the ref- 
erence system. We also stress the necessity to consider 
the chemical potential of the cluster as a variational pa- 
rameter. In Sec. IIIII we present numerical results for the 
phase diagram and the single-particle Green function in a 
one-dimensional (ID) system. Even in the simplest case 
of a reference system consisting of single-site clusters, 
the VCPT improves drastically on the CPT. In particu- 
lar, we obtain the correct form of the Mott lobes in the 
(t/U, /j,/U) phase diagram (t is the hopping amplitude, 
U the onsite repulsion, and /i the chemical potential). 
We find a good quantitative agreement with the results 
of the density-matrix renormalization groupAS when the 
number of sites in the cluster becomes of order 10. The 
last section is devoted to a summary of our results and 
a discussion of future developments, in particular the ex- 
tension of the VCPT to the superfluid phase. 
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II. VARIATIONAL CLUSTER PERTURBATION 
THEORY FOR BOSONS 

The Bose-Hubbard model is defined by the Hamilto- 
nian 

H = - ^(^t rr /^ r '+Hx.)-//^n r + — ^n r (n r -l), 

r,r' r r 

(1) 

where ip T ,ip\ are annihilation/creation bosonic operators 
and n r = ipj.ip T - The discrete variable r labels the sites of 
the lattice, which is assumed to be bipartite with co- 
ordination number z. The hopping matrix t satisfies 
t rr > = t > if r and r' are nearest neighbors and van- 
ishes otherwise; this assumption can however be easily re- 
laxed and longer-range hopping considered. For reasons 
explained below, the boson-boson interaction should be 
onsite. The density n, i.e. the average number of bosons 
per site, is fixed by the chemical potential [i. In the fol- 
lowing, we shall consider only the Mott phase and the 
zero-temperature limit. 



A. General formalism 

In the absence of superfluidity ((ip T ) = 0), the grand 
potential (per lattice site) fl and the single-particle Green 
function G can be obtained from the stationary point of 
the functional^fl»^i*^ 

n[G] = W l Trln (- G_1 ) + MGo'G - 1) + $[G]} , 

(2) 

where Go is the non-interacting Green function and $[G] 
the Luttinger-Ward functional. f3 — 1/T is the inverse 
temperature and N the number of lattice sites. Tr de- 
notes a trace over space and time indices. The stationary 
condition 5Q[G]/SG = yields the Dyson equation 

G = Gq E, 

where the self-energy E is defined by 

6$[G] 
13 ~ SG Tl • 



(3) 



(4) 



i and j are collective indices that label position and time. 
The SFA is based on the functional^ 

= ± {TrlnC-Go 1 + E) + F[E]} , (5) 



where 



F[E] = $[G] + Tr(EG) 



(6) 



is the Legendre transform of $[G]. In Eq. Jfjjl, G should 
be considered as a functional G[E] of the self-energy ob- 
tained by inverting QJ. F satisfies 



(7) 



and the stationary condition <5Q[E]/<5E = reproduces 
the Dyson equation J3J. 

So far, we have followed the approach of Ref. with 
minor modifications due to the fact that we consider 
bosons instead of fermions. We now introduce a refine- 
ment of the approach, the motivation of which will be 
discussed below. In a translation invariant system, the 
actual self-energy is diagonal in k space: E(k, k',z) = 
<5k,k'E(k, k, z) (z is a complex frequency). Without 
changing its value at the stationary point, we can there- 
fore modify the functional (0 into 



ft[E] 



1 

Wj3 



{Trln(-G 7 1 +E) + f[E]}, 



where 



E(k, k',z) = 6 k: k'E(k, k, z) 



(8) 



(9) 



is the diagonal part of the self-energy. One easily verifies 
that the stationary condition <5fi[E]/<5E = yields the 
Dyson equation 



G-\K z) = Go 1 (k,z)-E(k,z), 



(10) 



where G (k, iui v 



fi — ek- £k is the Fourier 



transform of —t r r / and gives the lattice dispersion of the 
bosons. In Eq. itlfill . E(k, z) denotes E(k, k, z), etc. Thus 
the two functional J^J and ® contain the same infor- 
mation in the case of translation invariant systems. 
Let us now consider the two Hamiltonians^ 2 



H(x) = H (x)+Hu, 
H(x') = H (x')+Hu. 



(11) 



H(x) is the Hamiltonian of the Bose-Hubbard model 
[Eq. QJ], and H(x') that of a reference system. Both 
Hamiltonians are defined on the same lattice and share 
the same interaction Hamiltonian Hjj . The kinetic parts 
Hq(x) and Hq(x') include a chemical potential term, x 
stands for the parameters on which Hq depends: the 
intersite hopping matrix t and the chemical potential 
(j,. The Hamiltonian H{x) is translation invariant, but 
that of the reference system, H(x'), may not be. The 
Luttinger-Ward functional $[G] is given by the sum 
of the two-particle irreducible (skeleton) diagrams^A*^ 
and is independent of Go- It follows that H{x) and 
H(x') share the same Luttinger-Ward functional $[G] 
and therefore the same functional F[E}22. This leads us 
to consider the functionals 

ttx[E] = J^TrlnHV+EHFp]}, 

MS] = ^{Trln(-G - 1 + E)+F[E]}, (12) 

where G = G (x) and G' Q = G (x') are the Green func- 
tions corresponding to H (x) and H (x'), respectively. 
Note that we have taken advantage of the translation 



3 



invariance of H(x) in writing O^fS]. The unknown func- 
tional F[E] can now be eliminated by taking the differ- 
ence of the two preceeding equations, 

ryq = rvp] + ^_{Trin(-Go 1 + 1) 

-Trln^^^ + S)}. (13) 

The form 111 3t of the self-energy functional f^E] is 
exact. To make the determination of the stationary 
points of fi x [E] possible, one has to restrict the space 
of self-energies. A natural approximation is to evalu- 
ate fi x [E] for the (physical) self-energies E(a;') obtained 
from the reference system - assuming that the refer- 
ence Hamiltonian can be solved exactly>22, The functional 
Q x [E(x')] = Q x (x') becomes a function of x', 

n x (x>) = n' + ^{TrlnhGo- 1 + 

- Trln[-Go _1 + E(x')]} (14) 

where O' = n x '[E(a;')] is the exact grand potential of the 
reference system. The stationary condition becomes 

dnjx') „ , N 

Note that the stationary point does not have to be a min- 
imum. We shall see that it is actually never a minimum 
when the chemical potential of the reference system is 
taken as a variational parameter. This conclusion also 
holds in fermion systems^ 

B. Reference system 

In the VCPT one considers a reference system consist- 
ing of a superlattice of N c = N/ L identical clusters with 
no intercluster kinetic coupling (L is the number of sites 
in a cluster) . The determination of fl' and the self-energy 
E(x') then requires to solve a finite size system Hamilto- 
nian provided that the interaction is local. For moderate 
values of L, this can be done numerically. Eq. 11411 can 
be rewritten as 

n x (x') = tf- J- £ ln[-G(k,iw n )]e <w »* 

k,u; n 

+ ^E trln [-G'^„)] e ^" (16) 

(u) n is a bosonic Matsubara frequency), where G _1 = 
Gq 1 — E and G /_1 = Gq _1 — E (from now on, we do not 
write explicitely the x' dependence of Y,(x')). tr denotes 
the trace over space indices, tr ln(— G') can be evaluated 
by considering a single cluster and multiplying the result 
by N c = N/L to take into account the total number of 
clusters. When the chemical potentials [i and p! differ, 



the usual factor e luJnTI (ry — > + ) is necessary for the sum 
over uj n to converge. 

Let us now discuss the motivation for using the func- 
tional (jHJ rather than ©. The reference system has 
the periodicity of the cluster superlattice, but not that 
of the translation invariant system of interest. Except 
for single-site clusters, the self-energy E is not diagonal 
in k space. By making use of the functional iJBJ, we 
ensure the translation invariance of the Green function 
defined by Eq. IIOII . In applications of the VCPT to 
fermionic systems, the translation invariant Green func- 
tion of the system was approximated by the diagonal part 
of G = (Gq 1 — E) -1 . However, because E itself is not 
translation invariant, the resulting excitation spectrum 
exhibits spurious gaps that reflect the periodicity of the 
reference system. To some extent, this difficulty can be 
eliminated by introducing an artificial broadening of the 
energy states; technically this is achieved by choosing a 
sufficiently large value of r\. In bosonic systems, even 
when one allows for an energy broadening, these gaps 
are so pronounced that the excitation spectrum does not 
bear much physical meaning. 

The translation invariant self-energy E(k, z) turns out 
to be simply related to the cluster self-energy when pe- 
riodic boundary conditions - as allowed in the VCPT, 
since this simply assumes a particular choice of the (in 
principle variational) intersite hopping matrix r - are 
chosen for the cluster. Let us write the self-energy of the 
reference system as 

E(R + r a , R' + r b , z) = 5 RiR ,X ab {z), (17) 

where T, a i,(z) is the cluster self-energy obtained by solv- 
ing numerically the Hamiltonian H(x'). We use the no- 
tation R + r a for a site of a lattice, R for a site of the 
superlattice, and r a (a e [1,^]) for the position within 
the cluster. The translation invariant self-energy then 
reads 

1 L 

E(k,z) = z e- 4k - (r "- rb) £ afc (z). (18) 

a,6=l 

With periodic boundary conditions for the cluster, we 
can define the Fourier transform of E a fc(z) for any vector 
k of the reciprocal superlattice. In that case, E(k, z) 
defined by Eq. ltl8|l is nothing but the interpolation of 
the translation invariant cluster self-energy to all vectors 
k of the Brillouin zone. 



C. Numerical implementation 

To numerically evaluate the sum over Matsubara fre- 
quencies in Eq. ltl6|l . one has to get rid of the conver- 
gence factor e luJnri . To this end, one subtracts and adds 
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the term 



5> cln 



Ah (K) 



-/3(y+A«(K)) 



lnfl 



(19) 



where tr c denotes the trace over intracluster space in- 
dices. Ahg = ho — h' is defined from the Green functions 
Gq 1 (iu> n ) = iuj n -h and Go _1 (iw„) = iu) n -h' Q . Ah (K) 
is the Fourier transform of Aho with respect to the su- 
perlattice of clusters, 

A/i (K) a6 =^e- <K -( R - R ')A^(R + r O) R' + r i ,). (20) 

R 

A/io(r,r') denotes the matrix elements of Aho in real 
space, and K is a vector of the reduced Brillouin zone 
corresponding to the superlattice. A Q (K) (a S [l,L]) 
are the eigenvalues of Aft-o(K), and the real number y 
should satisfy y > max a: K |A a (K)|, Eq. II19JI is derived 
in Appendix El Combining 11 fit and l|19ll , we finally 
obtain 

n x (x') = - ^E{E ln t- G ( k '^™)] 



tr ln[— G'(iu> n )] + tr c In 

K 

,-/%+A a (K)) 



1 - 



Ah (K) 



iu n - y 



-0V\ 



(21) 



The term inside the curly brakets behaves as l/w„ when 
\uj n \ — > oo, so that the convergence factor e luJ " v is not 
necessary anymore. Eq. l|2"T]) is the starting point of the 
numerical calculations discussed in Sec. IIIII 



D. /i' as a variational parameter 

In this section, we discuss the role of the chemical po- 
tential p! of the cluster as a variational parameter. The 
boson density is obtained from 



dn x (x') 



(22) 



where the value of x' — x'(x) is determined from the 
stationary condition 11 51 . Combining Eqs. 1221) and (11511 . 
we obtain 



dn x (x') 



dn x { x ') 

dx' 



dxf_ 



dn x {x') 



(23) 



(24) 



From Eq. llfill . we then deduce 

Here we have used the fact that the only dependence on 
/j, (at fixed x') comes from Go(k, iuj n ). Thus the station- 
ary condition (1 1 5 1 ensures that the approach is thermo- 
dynamically consistent: the boson density can be calcu- 
lated either from the grand potential or the single-particle 
Green function. 

The stationary condition also ensures that the boson 
density n of the system is the same as that of the reference 
system (n c ). To see this, we make use of the following 
result derived in Appendix IB1 



-/3Z T (k) 



^ln[-G(k,^„)K"" 

k,u> n 

-^lnll-e-^(k) 
^trln[-G , («w„)]e i "" ?? = 

- In 1 1 - ^ K -< I + E ln 1 1 - e ~ PKn 



EH 1 - 

k,7 



(25) 



a, 7 



where E 1 {\&) and Z 7 (k) [E' aj and Z'] denote the poles 
and the zeros of the Green function G(z) [G'(z)]. By 
virtue of the definition of the self-energy, Z 7 (k) and Z' 

correspond to the poles of the £( z ) and T,(z), respec- 
tively. We show in Appendix [0 that £(z) and E(z) 
share the same poles, i.e. {Z 7 (k)} = {Z' aj }. From 
Eqs. 111612511 . we then deduce 



k,7 



(26) 



a, 7 



(9 is the Heaviside step function) in the zero temperature 
limit. Note that the energies E' ai are N c times degener- 
ate, since G' is the Green function of the whole reference 
system (N c clusters). Eq. (l26l) gives the boson density 



N ^ 



dE 7 (k) 



k , 7 



-s 7 (k) 



(27) 



To relate n to the boson density n c = —dfl'/dp' in 
the reference system, we use the stationary condition 

dn x {x')/dii' = o, 



— V 

N ^ 



<9£ 7 (k) 



k,7 



dfx' 



i x ^ dE' a . 

a, 7 



?[-£? 7 (k)] 
9[-E' ay ]. (28) 
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This equation can be simplified by noting that the chem- 
ical potential fi' is a mere shift of the excitation energies 
E' so that E' + fi' is independent of //, 



dE. 



ay 



dfi' 



= -1. 



(29) 



Consider now the relation (1 1 Oil between G(z), Gq(z) and 
In this equation, E(z) is a function of z + fi', while 
Go(z) and G(z) are functions of z + \i. It follows that 
£" 7 (k) + n is a function of /i — //, which gives 



d£ 7 (k) 



<9£: 7 (k) 



From Eqs. p7!3()| . we deduce 



n = n c + 



a, 7 



(30) 



(31) 



k,7 



In the Mott phase, the coupling between clusters will 
transform the discrete excitation energies E' aj into en- 
ergy bands £7 7 (k). Since the gap in the excitation spec- 
trum remains nonzero, the number of negative energies 
cannot change, so that n = n c . 

Clearly, and this is confirmed by our numerical calcu- 
lations, the important point here is to take the chemical 
potential y! of the cluster as a variational parameter. 
Other parameters, such as the intersite hopping matrix 
t', do not have to be considered as variational in order 
to obtain a consistent result for the boson density. This 
point has also been stressed in Ref. HE 



III. NUMERICAL RESULTS IN ONE 
DIMENSION 

In this section, we present numerical results for ID 
systems. 

The simplest reference system consists of single-site 
clusters (L = 1). The only variational parameter is 
then the chemical potential x' = //. A single-site clus- 
ter can be solved analytically exactly. The state with 
p > particles is an eigenstate with the energy e p — 
—li'p + (U /2)p(p — 1). This yields the partition function 



-/30' 



p=0 



a. 
+ 

a 
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FIG. 1: Grand potential Qx(ji') vs. A// = // — fl in ID for 
t = 0.1 (single-site clusters). In all figures U = 1. 



The self-energy E(iw n ) = iui n + fx' — G' 1 (iu) n ) is local 
in real space and diagonal in reciprocal space: E = S, 
From Eq. IjlOJl and G^ik, ioj n ) = iu) n + fj,— ek, we obtain 



G(k,iui n ) = 



G'{iu n ) 



1 - (A// + e k )G'(iu n ) 
1 - Zk z k 



luJ, 



lU!. n 



(34) 



where 



Et 



-fi + - (-A/i' + e k + U) 



± - [U 2 + 6U(A[i' + Cfc ) + (Ap' + e k f] 



21 1/2 



(35) 



and Ap' = // — fi. Given that —fi',E k < and U — 

J k 



n',E£ > (see below), Eq. gives 



n x (p') = -2ii'-±J2 E k- 



(36) 



-0e p 



(T — > 0). (32) of the Mott gar 



For A/i' = 0, we recover the single-site CPT results 
obtained earlier i 1 ^ 22 ! 2 ^ A finite A/i' does not change the 
structure of the Green function G(k, iu n ). The excitation 
spectrum reveals the generic characteristics of the MI-SF 
transition^ There are two excitation branches, E^ > 
and E^ < 0, which coincide with the cluster excitation 
energies —fi' and U — fi' in the limit t — > 0. The dis- 
persion of Ej^ increases with t, which leads to a decrease 



Eto 



E 



k=0- 



(See, for instance, the 



The number of bosons n c in the ground state is obtained 
from e„ c = min p e p . This condition leads to n c — 1 < 
fi'/U < n c if fi' > —U, and n c = if fi' < 0. In the 
following, we consider only the MI with one boson per 
site, n = n c = 1, which requires < fi' < U. The zero- 
temperature Green function is local in space and reads 23 



G'tiuin 



7 + - 



[i'-U 



(33) 



figures in Ref. [2j for the case Afi' = 0.) By varying the 
chemical potential fi (with t and U fixed), one induces 
a transition from the commensurate incompressible MI 
to the incommensurate compressible SF when E^_ or 
E k= o vanishes. This transition is density driven and its 
critical behavior mean-field like. At the tip of the Mott 
lobe shown in Fig. E^_ and E^ =0 vanish simultane- 
ously, i.e. the Mott gap closes. This transition occurs at 
fixed density; it is driven by phase fluctuations and is in 
the universality class of the (d+ 1) XY model (with d the 
space dimension) i 
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0.2 0.4 0.6 0.8 1 

t 



FIG. 2: Phase diagram in ID obtained from a single-site clus- 
ter. The solid lines show the boundaries between the n = 1 
Mott insulator (MI) and the surrounding superfluid phase. 
The dashed line is the CPT result. The dotted lines show the 
values of p! at the Mott-insulator-to-superfluid transition. 
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FIG. 3: (color online) Phase diagram in ID obtained from a 
L-site cluster. The (red) squares show the results from the 
density-matrix renormalization group^ft 



Within the VCPT, one has to determine Ap' from the 
stationary condition dfl x (p')/dp' = 0, 



W + Ap' + e k 



N 



t [U 2 + 6U(Ap> + e*) + {Ap' + e k ) 



2ll/2 



3. (37) 



Using Eq. ffiTjl . one can verify that the boson density 
n = —dVt/dp obtained from Eq. l.ifill equals that of 
the cluster, n = n c = 1, in agreement with the gen- 
eral proof given in Sec. Ill PI In order for the square 
root in Eqs. H35I37I) to be defined, A// should satisfy 
A// < D - (3 + 2V2)U or A/i' > D - (3 - 2y/2)U, where 
D = — £fe=o = zt. (Note that the Mott gap vanishes when 
A// = D - (3 ± 2y/2)U.) The grand potential Q x (//) is 
shown in Fig.^for t/U — 0.1. We see that the stationary 
point of Q x (p') is a maximum with respect to variations 
of p! . This turns out to be always true regardless of the 
space dimension or the number of sites in the cluster. 
Since we expect n x (x') to be a minimum with respect 
to other variational parameters - for instance the intra- 
cluster hopping amplitude t' - the stationary point of the 
grand potential Q x (x') will be in general a saddle point 
for larger clusters (L > 2). 

Once A/i' is known, the values of the chemical poten- 
tial /i for which the Mott state becomes unstable are 
obtained from 2^L = and E^ =0 — 0. The low- 
est value of A/^' shown in Fig. corresponds to A/i' = 
D-(3-2V2)U. Ast/U increases, this value moves closer 
to the stationary point, and coincides with it when the 
Mott gap closes. This signals a transition to the super- 
fluid state with a density n = 1. For larger values of t, 
the ground state is superfluid for any value of /i (provided 
the boson density remains finite). 

The phase diagram of the ID Bose-Hubbard model in 
Fig- shows that the results obtained from the CPT and 
the VCPT differ drastically even for single-site clusters. 
Whereas the CPT gives the round-shape Mott lobe char- 
acteristic of mean-field theories - and independent of the 
dimension but for a trivial dependence on the number 
z of nearest neighbors -, the VCPT qualitatively repro- 



duces the shape of the Mott lobe in lDi Vfti 2 ^ In partic- 
ular, it yields a reentrant behavior of the MI state as t 
is increased with /i ~ 0.05J7 fixed, and a very pointed 
shape around the lobe tip. For sufficiently large values 
of t, p' comes near U/2, while the decrease of /x appears 
to be tied to the bottom of the free boson dispersion e^, 
i.e. /i ~ —2t + const. The gap closes very slowly as t 
increases, and the MI disappears for t/U ~ 2 (not shown 
in Fig.|2j). 

Note that the pointed shape of the Mott lobe in ID is 
usually attributed to the slow decrease of the Mott gap, 



E, 



k=0 



E, 



fe=0 



exp 



const 



t 



(38) 



near the Berezinskii-Kosterlitz-Thouless (BKT) transi- 
tion taking place at the lobe tip (t = t c ) i 10 i 26 The fact 
that the single-site VCPT reproduces the correct over- 
all shape of the Mott tip, while the physics of the BKT 
transition is clearly out of reach of this method, is quite 
remarkable. 

For larger clusters (L > 2), the ground state and the 
grand potential Q', as well as the single-particle Green 
function G", can be obtained numerically using the Lanc- 
zos method^ The intracluster hopping amplitude t' = t 
is held fixed, and only the chemical potential p! is taken 
as a variational parameter. As in the case L = 1, the 
stationary point dQ, x {n')/ 'dp! — is a maximum. For a 
given /x, not all values of p! are physically acceptable. On 
the one hand, p! has to be such that n c = 1 (we restrict 
ourselves to the n — 1 MI). On the other hand, we have 
to ensure that the ground state is stable. For L = 1, we 
saw that for certain values of p! , the single-particle exci- 
tation energies are complex. This instability also shows 
up as a disagreement between the boson density n ob- 
tained from the trace of G [Eq. jSUl ] and n c . We use this 
latter criterion to verify the stability of the ground state, 
since it does not require to obtain the excitation energies 
while minimizing the grand potential with respect to p' . 

The phase diagram of the ID Bose-Hubbard model is 
shown in Fig. 03 together with the results obtained from 
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FIG. 4: (color online) Spectral function A(k, uf) — 
Im.G(k,u)+ir)) obtained from 4-site clusters: t/U = 0.05, 
0.1, 0.2 and 0.3 (from top to bottom). 



the density-matrix renormalization group (DMRG). 1 ^ 
Single-site clusters give a good approximation to the up- 
per boundary of the n — 1 MI. By contrast, the lower 
boundary obtained within the VCPT strongly depends 
on the number of sites in the cluster. For L > 8, the 
VCPT reproduces fairly accurately the results obtained 
from the DMRG. 

We show in Fig. 31 the spectral function A(k,oj) — 
— 7r ImG(fc, cj + irj) (uu real). The spectrum consists of 
two well defined dispersion branches separated by the 
Mott gap. We see clearly how the Mott gap closes as we 
move closer to the Mott lobe tip. Note that besides the 
two main excitation energies, the spectral function ex- 
hibits additional (weaker) structures at positive energies 
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FIG. 5: (color online) Comparison between the single-particle 
density of states (u) > 0) obtained from the translation invari- 
ant self-energy E(fc, z) (upper (green) curve p(u>)) with that 
obtained from the self-energy T,(k,k',z) (lower (red) curve 
—p{ui)). The plots are obtained from 4-site clusters. 



when the Mott gap is small (see the two bottom graphs 
in Fig. QJ ; we do not know whether these have a true 
physical meaning or are due to the finite accuracy of our 
numerical calculations. 

As discussed in Sec. [H] by making use of the transla- 
tion invariant self-energy S, we obtain a boson dispersion 
free of spurious gaps coming from the periodicity of the 
reference system (Fig. 0J. This is illustrated in Fig. 
where we compare the single-particle density of states 
p(uj) = J 7j^A(k,uj) with that obtained from the non- 
translation-invariant self-energy S(fc, fc', z). In the latter 
case, two gaps arising from the periodicity of the 4-site 
clusters are clearly visible. 



IV. CONCLUSION 

The VCPT )^ 2 '^'^ 4 '^ 1 ? which was previously applied to 
strongly-correlated fermion systems ^2ML2M1 can be 
extended to boson systems. We propose a modification 
of the original formulation which ensures that the fi- 
nal results are translation invariant despite the fact that 
the reference system breaks translation invariance. This 
translation invariant VCPT is not restricted to boson sys- 
tems but should also apply to fermion systems. 

The results obtained for the ID Bose-Hubbard model 
indicate that the VCPT is an efficient method for study- 
ing the Mott-insulator-to-superfluid transition in boson 
systems. We stress the importance of taking the chemi- 
cal potential p! of the cluster as a variational parameter. 
This ensures that the approach is thermodynamically 
consistent (n = -dQ/dp, = -Tr(G)/(N0)), and that the 
boson density n is the same in the system and the refer- 
ence system. The grand potential is found to be a maxi- 
mum with respect to variation of //, which implies that in 
general it will correspond to a saddle point when several 
variational parameters are considered. Even for one-site 
clusters, the VCPT and CPT - where //' = (X is not a 
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variational parameter - differ drastically. Whereas the 
CPT gives the usual dimension-independent mean-field 
results, the VCPT reproduces the characteristic pointed 
shape and reentrant behavior of the Mott lobe in ID. 

The extension of our results to higher dimension does 
not raise any difficulty and will be discussed elsewhere. 
More subtle is the application of the VCPT to the su- 
perfluid phase. In order to describe a broken gauge- 
symmetry phase ((Vv) ^ 0), one should add to the ref- 
erence system Hamiltonian a field h' r that couples to the 
boson operator ip r . The set x' of variational parameters 
will therefore at least include the chemical potential // 
and the field h! , A finite value of b! at the stationary 
point implies a nonzero condensate (ip r ) and superfluid- 
ity. Whether or not the VCPT will satisfy the Hugenholz- 
Pines theorem and thus correctly describe the gapless Bo- 
goliubov sound mode should determine the applicability 
of the VCPT to the superfluid phase. 



Acknowledgments 

W. K. thanks the EPSRC for support through the 
Grant GR/S18571/01. 




FIG. 6: Integration contour used in Eq. ijAlfl . The black dots 
denote the bosonic Matsubara frequencies iui n . 



integration by part gives 
de 1 



2inf3 



ln(l- e -^) 



E 

Q = l 

L 



d_ 



In 1 - 



Aq(K) 

e + irj — y 
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APPENDIX A: DERIVATION OF EQ. JT9j| 



Eq. O follows from JUJ and ijAl)) . 



Let us consider the integral 



j^n B {z)tr c \n (l 



A^o(K) 



dz 



Aq(K) 

z-y 



(Al) 



where A Q (K) (a & [1, L]) are the eigenvalues of A/io(K). 
n B (z) = (eP z — 1) _1 is the Bose-Einstein distribution 
function. If we choose y > max Q: K |A Q (K)|, the branch 
cut of the logarithm in Eq. l|Alll does not extend to the 
origin z — 0, and we can evaluate I using the contour 
shown in Fig. where < A < y — max aj K |A Q (K)|. 
From the residue theorem, we then obtain 



/ = 



L 



Ja 2*7r ^ 



In 1 - 



Aq(K) 

e + 17] — y 



c.c. 



(A2) 



APPENDIX B: DERIVATION OF EQS. l(25|l 

Let us consider the integral 

dz 



I = 



2iir 



n B (z)trln[-G'(z)]e z?) 



= <f^n B (z)J2H-G' a (z)]e^, (Bl) 



where G' a {z) (a G [1,-^]) are the eigenvalues of G'{z). 
The (exact) bosonic Green function can be written as a 
sum of simple poles, 



(B2) 



where sgn(r a7 ) = sgn(E' ay ). The excitation energies E' aj 
are nonzero in the Mott phase. As a result 



-G' a (z = 0)=J2w L>0 



(B3) 



The factor n B (z)e zv ensures that the part of the con- 
tour at infinity does not contribute to the integral. An 



and the branch cut of the logarithm in Eq. (|B1J) does not 
extend to the origin. We can therefore evaluate / using 
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FIG. 7: Integration contour used in Eq. iBlt . 



the contour shown in Fig.0 



1 = \ E H-G' a (iu 




2iir 



n B (e) 



J2{H-G' a (e + iv)]-cc} 



The second of Eqs. ffijl follows from JB4tl and (IB5I . 

The first one can be derived similarly. Analog results 
for fermionic systems can be found in Ref. 133 . 



APPENDIX C: CLUSTER SELF-ENERGY E(z) 

In this appendix, we discuss two important properties 
of the cluster self-energy T,(z). 

1. lim^i^oo E(z) 

The cluster self-energy £(z) tends to a constant value 
equal the Hartree-Fock value 2Un when \z\ — > oo. A 
similar result holds for the fermionic Hubbard model: see 
for instance Appendix A in Ref. 0; the derivation given 
in this paper can be straightforwardly extended to the 
Bose-Hubbard model. 



--hill-e^l 



X J 2ztt (3 

^[G^ 1 ^ + ir,)d t G' a {e + ir,) - c.c], (B4) 



with Z' ai the zeros 



where < A < mhi an {\E' a 

of G'{z). The last line in llB4t is obtained by integrat- 
ing by part. The factor inside the last brackets vanishes 
unless e is near a pole or a zero of G' a (e). Near a pole, 
G' a (e) ~ r aj /(e + ir]- E' ai ), and [•••] = 2inS(e - E' ai ), 
The zeros of G' a {z) are given by the poles of the self- 
energy E(z), which shares the same analytical properties 
as G'(z) and can therefore be written as a sum of simple 
poles as in Eq. JB2ll (see Appendix [0 for a more detailed 
discussion). It follows that G' a (z) has only simple zeros 

s a ~j(z Z 



and can be approximated by G' a (z) 
:tor inside the brai 
Z'), We deduce 



near z — Z'„^, The factor inside the brackets in Eq. JB4I) 



a) ■ 

then gives —2iirS(e 



-iVlnll 
(3^ 



a, 7 



iVlnll-, 
(3^ I 



a, 7 



(B5) 



2. Poles of E(z) and E(z) 

The (matrix) self-energy £(z) shares the same analyt- 
ical properties as the Green function G'(z). Its Lehmann 
representation 

E(z) = E(oo) +£|a,7>-^_<a, 7 | (CI) 



(we use a bra-ket notation) leads to 



s(M) = (k|s(z)|k) 

= >:(oc)-£j<k. 



^Zi Q< ^ 



(C2) 



where S(oo) = 2Un. Since both {|k)} and {|a, 7)} span 
the whole Hilbert space, (k|a, 7) cannot vanish for all k. 
This implies that all poles Z ai of £(2) also show up as 
poles of S(z). 
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